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This paper presents a study of the two-stream instability of an electron beam propagating in a finite-size 
plasma placed between two electrodes. It is shown that the growth rate in such a system is much smaller 
than that of an infinite plasma or a finite size plasma with periodic boundary conditions. Even if the width of 
the plasma matches the resonance condition for a standing wave, a spatially growing wave is excited instead 
with the growth rate small compared to that of the standing wave in a periodic system. The approximate 
expression for this growth rate is 7 ~ {l/13)ujpe{nb/rip){Luipe/vb) ln(LLOpe/vb) [1 — O.lScos [Lujpe/vb + 7 >'/ 2 )], 
where ujpe is the electron plasma frequency, rib and rip are the beam and the plasma densities, respectively, Vb 
is the beam velocity, and L is the plasma width. The frequency, wave number and the spatial and temporal 
growth rates as functions of the plasma size exhibit band structure. The amplitude of saturation of the 
instability depends on the system length, not on the beam current. For short systems, the amplitude may 
exceed values predicted for infinite plasmas by more than an order of magnitude. 
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I. INTRODUCTION 


Interaction of electron beams with plasmas is of con¬ 
siderable importance for many plasma applications where 
electron emission occurs from surfaces. The electrons ex¬ 
tracted from the surface and accelerated in the sheath 
form a beam of electrons; the beam propagating in the 
plasma excites electron plasma waves through the two- 
stream instabilityji Laboratory plasmas and plasmas in 
industrial applications are usually bounded by electrodes. 
We show that electrodes greatly affect the growth of 
the two-stream instability. Though beam-plasma sys¬ 
tems have been studied extensively in the past using 
kinetic simulations,— the presence of electrically con¬ 
nected boundaries changes the character of the two- 
stream instability from convective to absolute, similar 
to the instability of a Pierce diode— In the Pierce diode, 
the instability was studied extensively taking only beam 
electrons and neutralizing ions into account as relevant to 
vacuum diodes, see e.g. Ref. 0 and the references within. 
Here, we consider the two-stream instability between a 
low density electron beam and high density plasma elec¬ 
trons as relevant to discharges. In this Letter, we have 
performed an analytical study and fluid and particle-in¬ 
cell simulations in order to obtain the growth rate of the 
two-stream instability in a finite plasma bounded by elec¬ 
trically connected electrodes. To the best of our knowl¬ 
edge and to some extent to our surprise the solution to 
this problem was not reported before. 

The linear stage of the instability can be described 
making use of fluid formalism which includes the conti¬ 
nuity equations 


drie^b dve,bne,b 

dt dx 


( 1 ) 


the momentum equations 


dVe,b , dve,b e 

nj. '^e,h o — -^5 

at OX m 


and the Poisson equation 
d^cj) 


dx‘^ 


= 47re (rie + rib — ni) 


( 2 ) 

( 3 ) 


where n^ b and Ve,b are the densities and the velocities 
of the plasma and beam electrons, —e and m are the 
electron charge and mass, E = —d(j>/dx is the electric 
field, (j) is the electric potential, and rii is the ion density. 
The initial plasma state is neutral: riefi + nbfi = riifi, 
where nf,,o and nb,o are the initial densities of the bulk 
and the beam electrons, and riifi is the initial density of 
ions, respectively. The ion density is uniform and con¬ 
stant, rii = rii^Q = const. Initially, the bulk and the 
beam electron densities and the beam flow velocity are 
uniform everywhere. Note that everywhere in this paper 
subscripts e and b denote values related to plasma and 
beam electrons, respectively. 

For the studies described in the present paper, the 
boundary conditions are non-periodic and describe a 
plasma produced in a discharge between two electrodes. 
At the ends of the system a; = 0 and x = L, the poten¬ 
tial perturbations are set to zero, (/>(0) = (j){L) — 0. The 
beam is injected at the boundary a: = 0. The bound¬ 
ary conditions for the beam electrons are nb{0) = Ub.o 
and Vb(0) = Vb,o, where Vb,o is the injection velocity of 
the beam. Note that in fluid simulations, a small sheath 
forms near the electrodes and more accurate boundary 
conditions are required to account for the sheath effect— 

The paper is organized as follows. In Section [TTl a dis¬ 
persion relation for the finite-length beam-plasma sys¬ 
tem is derived. In Section nni complex frequencies and 
wavenumbers obtained by direct solution of the disper¬ 
sion equation are compared with the fluid simulation 
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and approximate analytical formulas for the frequency, 
wavenumber, and temporal and spatial growth rates are 
given. Section Hvl compares growth rates in kinetic simu¬ 
lations with the predictions of the fluid theory. Section IVl 
discusses the amplitude of saturation of the instability 
and provides analytical formulas for the estimate of the 
saturation electric field amplitude. The results are sum¬ 
marized in Section [VII 


II. ANALYTICAL SOLUTION 

The dispersion equation is obtained by solving lin¬ 
earized Eqs. (PEI) for perturbations of plasma and beam 
electron densities and velocities. The perturbations are 
defined as 


dUg TLq ^e,0; ^b ^6,0: 

Svb = Vb- Vb^o, 5ve = Ve- 

Linearized equations can be readily solved using 
Laplace’s method.— However, we are only looking for 
an asymptotic solution which the system approaches on 
longer times. Following the Pierce methodj^ the asymp¬ 
totic solution for the potential has the following form: 

x) = {Ax + -f + D) 6 ““*, (4) 

where lo is the frequency of the wave, k± are the wave 
vectors of the two waves propagating in the system, and 
coefficients A, B, C, D are complex constants. The den¬ 
sity and the velocity perturbations are 

5ne,b{t,x) = (5n',^b + 

5ve,b{t,x) = (du;, + 6 ““*, 

( 5 ) 

The linearized equations for the parts of the perturba¬ 
tions proportional to exp(—iwt -I- ik±x) are 

—iojdnf + ik±dvfnf,Q = 0, 

—iujSv^ = — ik±6(j)^, 

m 

-iuSn^ + ik± [dvfubfi + VbfiSnf) = 0, 

{—iu} + ik±Vb o) = — ik±6(j)^, 

’ m 

—k'^Scj)'^ = iire {6nf + i5n^) , 
where Scj)'^ = B and S(j)~ = C. These equations yield: 


Substitution of relations m into the Poisson equation 
gives usual dispersion relation for waves 



(w 


LO 


2 

b,0 


- Vb,ok±)'^' 


( 7 ) 


Here = dTre^ng^g/m and w^g = 47 re^nb_o/™ are the 
electron plasma frequencies corresponding to the plasma 
and beam densities. 

The uniform parts of the density and velocity pertur¬ 
bations ®, which are proportional to exp(—itaf) and cor¬ 
respond to high-frequency uniform electric field given by 
the hrst term in Eq. o, are obtained in a similar way: 


? P 

Sv' = Sv{ = - A, 5n' = 5n[ = 0 . ( 8 ) 

ujm 

These perturbations correspond to high-frequency cur¬ 
rent flowing through the plasma and allow for 5ve ^ 0 
at the systems ends; Svb{0) = 0 because beam is injected 
with a given velocity but Svb{L) ^ 0 . 

Applying four boundary conditions Snb(O) = Svb(O) = 
S4>{0) = 54>{L) = 0 to perturbations (El and (0 and 
taking into account m and © in the form 


uj - k±Vb,o = ±' 


^b,0 


( 9 ) 


1 - 


^e ,0 


gives the following additional relation between uj and k: 


kl - 1 ) - 
kl - 1 ) 


ik'^k^wL 
oj - k+Vb,o 
ik'^k-UjL 

UJ - k-Vbfi' 


( 10 ) 


Eqs. ® and ([TOl determine the temporal [Im(a;)] and 
the spatial [Im(A;)] growth rates of the instability as well 
as the frequency [Re(a;)] and the wavenumber [Re(A:)]. If 
plasma electrons are absent and only beam electrons are 
taken into account (ne,o = 0), Eq. (flUl) reduces to the 
Pierce’s dispersion relation for vacuum diode. 

In order to solve the dispersion relation (flUl) , we intro¬ 
duce a new dimensionless variable 


X = 


^b,o/^e,0 



( 11 ) 


Substituting dm into m and assuming that uj = uJe,o in 
the left-hand side of (O gives 


k± 


(ITX) 


We ,0 

Vb,0 


( 12 ) 


Svf 


UJ 

k rie ,0 ^e ,0 

ca - Vb,ok± Sn^ 
k± ribfi 


- 


m uj^ 


nb,o 


(w - kv^) 

( 6 ) 


Substitution m into m yields equation for y 


-f - 1- 


(1 + x)x 
( 1 -x)^ ^ 
(1 + x)^ 


^l{l+X)Ln _ ^ 


(13) 


= 0 
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where L„ = LuJe,o/vbfi is the normalized gap width. 

Equation 0 gives X as a function of The fre¬ 
quency is calculated from m and for a low-density beam 
with rib^o <C ne,o it is 


U! = 


We,0 



^ 6,0 



(14) 


The wavenumbers k± can be obtained from 0. 

Function x(Tn) is complex with band structure, i.e. 
it changes abruptly at certain = c + 2 ttI, where c is 
a constant and I is an integer. Indeed, in the limit of 
Ln ^ 1, the first two terms in psp are dominant which 
gives the following approximate expression: 

-zxT„e-*x^" =2L2e-*^". 

Here, we used the fact that |x| "C 1 and Im(x) > 0. The 
solution of this equation is the Lambert or productlog 

functionjiS 


-ixL„ = W{2Lle-^^’'). (15) 

This function has many branches, the branch selected 
must ensure the maximal growth rate. When parameters 
of the plasma, e.g. the discharge gap, change, a transition 
from one branch to another may occur and the instability 
growth rate will change abruptly. 

Since x is complex and independent on it follows 
from m that the temporal growth rate of the instability 
is proportional to We,o(’T-b,o/’^e,o) unlike the growth rate 
of the resonant perturbation k ss cJe,o/vb,o in a periodic 
system proportional to uje,o (nb,o/'ne,o)^^^— ■ 

The analytical solution is verified by fluid and particle- 
in-cell (PIC) simulations described below. 


III. FLUID SIMULATIONS 

The fluid numerical model solves Eqs. ([IJ-ilS]). The 
densities in © are advanced using the SHASTA 
method. The velocities in ([5]) are advanced using an up¬ 
wind schemeji^ The model demonstrates excellent agree¬ 
ment with the theoryi in simulations of the instability 
of a cold beam in a cold plasma with periodic boundary 
conditions. 

The fluid simulations are carried out with the follow¬ 
ing common parameters: ne,o = 2 x 10^^ m“^, a;e,o = 
2.52 X 10^^ s“^, beam energy 50 eV and beam veloc¬ 
ity Vb,o = 4.2 X 10® m/s, the numerical grid cell size 
is 1.3 /im, the time step is 0.9 fs. The selected val¬ 
ues of spatial and temporal steps ensure stability of the 
SHASTA algorithm. The resonant beam wavelength 
Af, = 2TTVbfi/u!e,o is 1.044 mm for these plasma param¬ 
eters. Initially, the bulk electron flow velocity is given 
a harmonic perturbation Sve = 6ve,o siTi{xu!e,o/Vb,o) with 


E 

E 
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E 

E 


X 
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FIG. 1. Evolution of the bulk electron density perturba¬ 
tion in time and space in fluid simulation with L = lAt and 
a = 0.0006. Panels (a) and (b) correspond to the very begin¬ 
ning of the fluid simulation (a) and to the asymptotic constant 
growth stage (b); the corresponding temporal growth of the 
electric field amplitude is shown by the red curve in Fig. O 
Solid black lines in (a) and (b) represent propagation with the 
unperturbed beam velocity. Dashed black lines in (b) repre¬ 
sent phase velocity of the wave calculated as Re(a;)/Re(A:), 
where Re(a;) = 2.522 x 10^° s“^ and Re(fc) = 7.288 mm“^. 
Arrows A, B, and C mark times tA ~ 0.35 ns, ts ~ 3.01 ns, 
and tc = 141.8 ns when profiles shown in Figs. [2^, ^jp, and 
H; are obtained. 


the wavelength corresponding to the resonance in a peri¬ 
odic or an infinite plasma, the amplitude of the pertur¬ 
bation is very small, 6vefi = 0.1 m/s. 

The oscillations have the wavelength of the initial per¬ 
turbation during only the first few periods. The initial 
oscillation pattern corresponds to a standing wave. As 
the instability develops, the standing wave transforms to 
a propagating wave, see Fig. [T^. This process is accom¬ 
panied by the shrinking of the wavelength, compare the 
density perturbation profiles at three consecutive times 
in Fig. [2] At the initial phase of the instability, the per¬ 
turbations propagate with the original beam velocity, see 
Fig. [T^. At the asymptotic stage given by Eq.(|3]) with 
the spatial growth rate along the beam propagation, the 
wave phase velocity is noticeably lower than the veloc¬ 
ity of beam propagation, compare the slope of the black 
dashed line with that of the black solid lines in Fig. [TJd. 

Simulation reveals that before the asymptotic state es¬ 
tablishes, the temporal growth rate changes with time, 
see Fig. [21 Initially, the growth rate is large compared to 
the analytical value defined by Eqs. (1141) and (ITSl) . Then 
it gradually decreases towards the asymptotic value pre¬ 
dicted by the theory and it stays approximately constant 
for tens and even hundreds of plasma periods until the 
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FIG. 2. Spatial profiles of bulk electron density perturbation 
obtained at Ia ~ 0.35 ns (a), Ib ~ 3.01 ns (b), and tc = 
141.8 ns (c). Times tA,B,c are shown by arrows A, B, and C 
in Fig. [T] 



for 20 s < i < 160 s. In simulation sets one and two, for 
some values of L such a stage never appears, see the blue 
curve in Fig.[3l These values of L correspond to the gaps 
in the simulation data seen in Fig. 0] 

Overall, there is an excellent agreement between the 
simulations and the theory. The dimensionless values 
of [Re(w) - We,o]/(we,oa), Iin(w)/(we,oa), Re(fcA&), and 
Im(fcAb) obtained in simulation sets one and two (red and 
black curves in Fig. |4]) and by analytical solution of the 
theoretical dispersion relation (blue crosses in Fig.|T]) are 
very close to each other and appear to be functions of the 
dimensionless system length only, as predicted by the an¬ 
alytical solution given by Eqs. (USD and (HI- These func¬ 
tions for Re(a;), Im(w), and Re(fc) have band structure. 
Mathematically, it is the consequence of the presence of 
many branches in the Lambert function. The instability 
growth is given by the maximum growth rate value that 
changes from branch to branch when the gap size crosses 
some critical value, typically when L/Ab approaches an 
integer, see Fig. H) Similar band structure was also ob¬ 
served for the Pierce diodei^d^ FigureH^ shows the num¬ 
ber of wave periods in the gap as a function of the gap 
length. In all cases, it is very close to an integer number, 
although not exactly: 

Re(fc)L/(27r) ~ [L/Abl, (16) 

where \x\ = ceiling{x) is the smallest integer not less 
than X. 

Since the shape of the functions is universal for various 
beam densities, it is reasonable to introduce approximate 
formulas which fit the numerical solution as follows: 

Re(w) Ri ln(L„) [1 - 0.9 cos (L„ + 0.4)], (17) 

io 


FIG. 3. Amplitude of electric field oscillations versus time in 
fluid simulations with a = 0.0006 and L = 4A6 (red curve), 
L = 4.7A(, (blue curve). The curves are obtained in the point 
with coordinate x = 3.55 mm (red) and x = 4.34 mm (blue). 


nonlinear stage of instability and its saturation occurs, 
see the red curve in Fig. |3l Note that the modification of 
the wavelength mentioned above stops when the instabil¬ 
ity reaches the asymptotic stage, which for the red curve 
in Fig. [3] occurs near 20 ns. 

In order to investigate the dependence of the growth 
rate on plasma parameters, four simulation sets are dis¬ 
cussed below. In set one, the ratio of the beam to plasma 
density is a = nb^/up^o = 0.00015, the size of the sys¬ 
tem L increases from At, to 8.5At,. Set two is similar to 
set one but the beam density is higher, a = 0.0006. In 
set three, L — 3.4Ab is constant while a changes from 
0.0001 to 0.0006. The fourth set is similar to set three 
but L = 8.3Af,. 

In all simulations, the growth rates, the frequencies, 
and the wavenumbers are calculated during the asymp¬ 
totic stage when the temporal growth rate is constant for 
a prolonged period of time, see the red curve in Fig. [3] 


Im(a;) 


We.OO! 

13 


Tn ln(7/n) 


1 — 0.18 cos 



Re(fc) 


U1e,0 2 1 + 2.5cOs(£n) 

^6,0 _ 1-lT^ 


1 

(18) 

(19) 


In,(fc) « ^2HLn)-0.5_ ^20) 

The wavenumber and the spatial growth rate depend 
on the system length but are virtually insensitive to the 
beam density, see Fig. [5t and Fig. [5jl. The temporal 
growth rate is approximately linearly proportional to the 
relative beam density a. The linear law holds especially 
well for short systems, see the red curve in Fig. [5 ]d and 
compare red and black curves for L/Af, < 6 in Fig. IDd. 
For longer systems, however, deviation from the linear 
law becomes noticeable. 


IV. TEMPORAL GROWTH RATE IN KINETIC 
SIMULATION 

Kinetic simulations are carried out with the EDIPIC 
1D3V particle-in-cell (PIC) codeJ^ The code is modified 
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FIG. 5. Frequency (a), temporal growth rate (b), wavenumber 
(c), and spatial growth rate (d) versus the ratio of the beam 
and bulk electron densities in fluid simulations with L = 3.4At, 
(red) and L = 8.3Ai, (black). In (c), the black dashed line 
marks the resonant wavenumber. 


FIG. 4. Frequency (a), temporal growth rate (b), wavenum¬ 
ber (c), spatial growth rate (d), and the number of wave 
periods per system length (e) versus the length of the sys¬ 
tem. The blue crosses mark values obtained by analytical 
solution given by equations (HU, m, and ca. Solid red 
and black curves represent values obtained in fluid simula¬ 
tions with a = 0.00015 (red) and a = 0.0006 (black). Solid 
green curves are values provided by fitting formulas (HZl), (HU, 
HU}, and (1201) . In (c), the black dashed line marks the reso¬ 
nant wavenumber. 

to reproduce conditions of the fluid simulations. The 
ions form an immobile background, the boundaries have 
zero potential. The bulk electrons are reflected specu¬ 
larly from the boundaries. The beam electrons penetrate 
through the boundaries freely. The initial plasma den¬ 
sity and the beam energy are the same as in the fluid 
simulations. Collisions are omitted. Two simulations are 
carried out with L = 8.3Ah, a = 0.0006 but different 
number of particles per cell. One simulation has 10000 
particles per cell. The other simulation has 2000 parti¬ 
cles per cell. Below these simulations are referred to as 
10 k and 2k simulations, respectively. 

PIC simulations start with a significant level of statis¬ 
tical noise which is few orders of magnitude higher than 
the initial perturbation induced in the fluid simulations 


above. At the same time, the amplitudes of nonlinear 
saturation of the instability in PIC and fluid simulations 
are close to each other. As result, the time when the 
oscillations grow from the initial noise level to the satu¬ 
ration in PIC simulation is much shorter than that in 
a fluid simulation. Moreover, at the initial stage the 
growth rate gradually decreases which furthermore limits 
the duration of the asymptotic stage described by ana¬ 
lytic solution. For example, in the 10k simulation, the 
asymptotic stage lasts from 10 ns to 20 ns while in the 
fluid simulation that stage occurs between 10 ns and 45 
ns, compare the green and the red curves in Fig. [6] The 
short asymptotic stage in the 10k simulation still allows 
to calculate the temporal growth rate which appeared 
to be very close to the value obtained in fluid simula¬ 
tions. In the 2k simulation, however, the noise level is 
higher and the asymptotic stage is very short and barely 
detectable, see the blue curve in Fig. [SI 


V. SATURATION AMPLITUDE IN KINETIC 
SIMULATION 

PIC simulations described below are carried out with 
the following common parameters. The initial uniform 
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FIG. 6. Amplitude of electric field oscillations vs time in simu¬ 
lations with L — 8.3Xb and a = 0.0006. The curves represent 
fluid simulation (red), PIC simulation with 10000 particles 
per cell (green), and PIC simulation with 2000 particles per 
cell (blue). 


plasma electron density is rie.o = 2 x 10^^ m“^, the 
initial electron beam energy or the energy of injection 
is Wb = 50 eV, the beam-to-plasma density ratio is 
a = 1.5 X 10“^, the initial plasma electron temperature 
Tefi = 0.5 eV, the size of a cell of the computational grid 
is Ax = 2.078 X 10“® m corresponding to XD,e/8 where 
XD,e is the electron Debye length of the ambient plasma, 
both the plasma and the beam initially are represented 
by 2500 macroparticles per each cell of the grid. The 
ions are represented by an immobile uniform background 
with density which ensures that the plasma-beam system 
is initially neutral. The electron beam propagates in the 
positive x-direction. 

The following three PIC simulations are carried out. 
First simulation has periodic boundary conditions and 
the system length of L = 2Af, where Af, = 27rt4/we_o 
is the wavelength of the plasma wave resonant with the 
beam in an infinite plasma, 14 is the beam velocity. For 
the selected parameters, Xt = 1.044 mm. Second sim¬ 
ulation has non-periodic boundary conditions similar to 
the ones used in the fluid model. The boundaries are 
grounded, the plasma electrons are reflected from the 
boundaries while the beam electrons penetrate through 
them freely. System length in the second simulation is 
L = 8.3Af, which corresponds to the maximum of the tem¬ 
poral growth rate in the dispersion band with L/Xb ~ 8 in 
the fluid simulation. Third simulation also has the non¬ 
periodic similar to the second simulation, but the system 
length is much shorter, only L — 1.5A{,. In the second 
and the third simulations, the beam injection occurs at 
the boundary x = 0. 

The theory of interaction of a cold beam with a cold 
plasma predicts that the exponential growth of the am¬ 
plitude of plasma oscillations is followed by saturation 
and subsequent amplitude oscillationsSuch a picture 
is reproduced in the first simulation, see Figs. (iKa) and 
(c). The plasma wave propagates in the direction of beam 
propagation and has constant amplitude along the sys¬ 
tem, see Fig. [TKb) and Fig. EKa). The theoretical growth 
rate is 

( 21 ) 


and the electric field amplitude in the first maximum is 

^i,ma. = 3kWba^^^, (22) 

where We,o is the Langmuir frequency of plasma electrons 
and k = uJefilVb is the resonance wavenumber. For the 
selected beam and plasma parameters, a;e,o = 2.52 x 
10^° s“^, Vb = 4.19 X 10® m/s, and k = 6015.9 m“^. 
Therefore, the theoretical growth rate (I^Tl) is 

Im{uj) = 0.938 X 10® s“i 
and the electric field amplitude maximum (1221) is 
El,max = 2549V/m . 

Both values are very close to the simulation results, see 
the red curve in Fig.[7Kc) and compare it with the dashed 
straight line which corresponds to the theoretical growth 
rate ([?T|l . It is necessary to mention here that during 
the first 4 ns, the growing oscillations are obscured by 
the noise present in the system due to the finite number 
of particles in simulation. The saturation begins at t=6 
ns when the beam particles start passing each other, see 
the phase plane in Fig. EKb). Note that by this time the 
beam electrons in the laboratory frame travel about 25 
mm which is several times more than the length of the 
system in the second simulation. 

In the second simulation, the boundary conditions are 
non-periodic and the linear stage of the instability follows 
the fluid theory for the finite-length systems developed 
above - the wave amplitude grows both along the sys¬ 
tem and in time, the temporal growth rate is close to the 
fluid value, compare the red curve with the straight black 
dashed line in Fig.[7Kg). The saturation of the amplitude 
occurs around 70 ns and here the amplitude is maximal 
near the exit end of the system, see Fig. [T^d) and (f) 
near arrow 1 and the electric field profile in Fig. He). 
The maximum wave amplitude is an order of magnitude 
higher than that in the periodic simulation, and it causes 
much stronger velocity perturbations of the beam parti¬ 
cles. Due to both limited distance of interaction between 
the beam and the wave and the nonuniform wave ampli¬ 
tude, the beam electrons start passing each other only 
near the exit end, see Fig. EKd). 

An interesting process occurs after the first saturation. 
Position of the maximum amplitude gradually moves to¬ 
wards the injection boundary until it reaches the distance 
of about 5.5 mm, and the value of the maximum almost 
doubles, see Fig. [TKd) and Fig. EKe). After 150 ns, it is 
this new maximum where the passing of beam electrons 
is achieved, not the exit end of the system, see Fig. Elf)- 
Downstream of this maximum (x > 6 mm) the beam 
electrons are completely mixed in the phase plane. As 
a result of this change in the beam structure, while up¬ 
stream of the maximum (x < 5 mm) the wave propagates 
along the beam direction with spatially growing ampli¬ 
tude, downstream of the maximum the wave pattern is 
closer to that of a standing wave, compare Figs. [TKe) and 

(f)- 
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FIG. 7. Results of PIC simulation with periodic boundaries (a,b,c), non-periodic boundaries and L = 8.3A6 (d,e,f,g), and 
non-periodic boundaries and L = l.SAb (h,i,j). The top row (a,d,h) shows amplitude of oscillations versus coordinate and 
time. The middle row (b,e,f,i) shows electric field versus coordinate and time. The bottom row (c,g,j) shows the amplitude of 
oscillations versus time at certain locations marked by vertical arrows in (a), (d), and (h), respectively. The green and the red 
curves in (g) correspond to locations marked by vertical arrows A and B in (d), respectively. The horizontal arrows in (a,b) 
mark time of the snapshot shown in Figs. [8l^a,b). The horizontal arrows in (d,e,f) mark time of snapshots shown in Figs. |8l^c,d) 
(arrow 1) and in Figs. |8Ke,f) (arrow 2). The horizontal arrows in (h,i) mark time of the snapshot shown in Figs. [HKgjh). The 
dashed black straight line in (c) shows the exponential growth with the theoretical growth rate in an infinite plasma. The 
dashed black straight lines in (g) and (j) corresponds to the growth rates obtained in fluid simulations with the same system 
length and beam current. 


The second simulation clearly shows that the ampli¬ 
tude of saturation of the two stream instability in finite 
length plasmas can be significantly higher than that in an 
infinite plasma. In a bounded system, the length of in¬ 
teraction between a beam electron and the plasma wave 
cannot exceed the distance between the boundaries. An 
infinite plasma has no such a limit. Here a wave of mod¬ 
est intensity can interact with beam electrons over longer 
distances before the mixing of the beam electrons in the 
phase plane occurs. In order to achieve such a mixing 
over much shorter distances, which is the case in bounded 
systems, the wave field must be much stronger. 

In order to find how strong this effect can be, the third 
simulation with the non-periodic boundary conditions is 
performed with a very short system length L = 1.5A;,. 
The wave pattern corresponds to a standing wave with 5 
nodes and 4 antinodes. The antinodes and the 3 middle 
nodes are clearly visible in Figs. [Tjh) and (i). There are 
two nodes with approximately zero electric field at the 
ends of the system. These two nodes were not resolved 
by the code diagnostics used to produce Figs. [TKh) and 
(i), but they are visible in the electric field profile in 
Fig.[5Kg). The temporal growth rate in this simulation is 
about 34% higher than the growth rate in the fluid simu¬ 


lation with the same parameters, compare the red curve 
with the black straight dashed line in Fig.[3(j). The max¬ 
imal wave amplitude reaches 170 V/mm which is more 
than 60 times stronger than the field in the periodic sys¬ 
tem, compare Figs. [T^j) and (c). Such a strong electric 
field produces mixing of beam electrons on a very short 
distance of 1 mm which is about one resonance wave¬ 
length, see Fig.IHKh). 

The phase plots shown in Fig. [8] prove that in the fi¬ 
nite length system the saturation of the instability oc¬ 
curs when the beam particles are overtaking each other. 
This process depends on the wave amplitude and the 
system length but should not depend on the beam cur¬ 
rent. To check this, two additional simulations are car¬ 
ried out with L = 8.3Ab and the relative beam density 
of a = 0.0003 and 0.0006. Another two additional sim¬ 
ulations with these beam densities are carried out for 
L = 1.5Aft. The results of these simulations combined 
with those obtained above for a = 0.00015 confirm that 
the amplitude of the first maximum of saturation of the 
instability is virtually insensitive to the beam current, see 
the red and the green curves with markers in Fig. [Hlja). 
The only difference is that for higher current the satura¬ 
tion is achieved faster. Note that in the whole range of 
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FIG. 8. Electric field profiles (a,c,e,f) and electron “velocity versus coordinate” phase planes (b,d,f,h) in the PIC simulations 
with periodic boundaries (a,b), non-periodic boundaries and L = 8.3A(, (c,d,e,f), and non-periodic boundaries and L = 1.5Af, 
(g,h). Snapshots (a,b) are obtained at time 6 ns shown by the horizontal arrows in Figs. [7ta,b). Snapshots (c,d) are at time 
71.88 ns shown by arrows 1 in Figs. [T^d,!). Snapshots (e,f) are at time 185.69 ns shown by arrows 2 in Figs. [7td,e). Snapshots 
(g,h) are at time 979.8 ns shown by the horizontal arrows in Figs. [7th, i). In the phase planes, the beam electrons are represented 
by the blue color while the plasma electrons are represented by the color map, the white background correspond to the empty 
space. 


the beam density a considered, the values of the satura¬ 
tion amplitude in the finite system are much higher than 
the predictions for the infinite system, compare the red 
and the green curves with markers with the blue curve 
in Fig. HJa). 

The wave amplitude which causes overtaking of beam 
particles can be estimated as 


eE 

TOe(w — kVhY 


Ab 


(23) 


where the left-hand side is the displacement of the parti¬ 
cles trapped by the wave in the wave frame. Using equa¬ 
tion (EH) one can replace — with (ywe^o)^- Then, 
replacing XbUie,o with Vbfi one can write an expression for 
the maximal electric field of the wave as 


E 


max 


meVb^0^e,0 

e 


Ix^l- 


(24) 


For estimates, the value of is convenient to find as 


= •( 1 — -^^^[i?e(A:) -I- ilm{k)] 

We,0 


(25) 


with Re(k) and Im(k) given by the approximate formu¬ 
las EH) and EH. Note that expressions EH and EH 
are independent on beam current and are functions of 
the normalized plasma gap width only. Therefore, 
the maximal field (l24ll depends on the beam velocity Vb 
and the gap width L but does not depend on the beam 
current. 


A dependence Emax{L) calculated with (l25|l for the 
beam parameters used in the simulations above is shown 
by the black curve in Fig. IHb)- The oscillations in 
this curve reflect the band structure of the wave num¬ 
ber in the finite length system. The saturation values 
obtained in the PIC simulations are remarkably close to 
the values given by Eq. (|25|) , compare curves with mark¬ 
ers with the horizontal dashed curves of the same color in 
Fig.lHa), also compare the markers with the black curve 
in Fig. IHb). The value of E^ax decreases with L and 
eventually approaches the saturation values for the infi¬ 
nite system given by Eq. EH, compare the black curve 
with the horizontal blue lines in Eig. IHb)- 

PIC simulations discussed above demonstrate that the 
growth of the maximal electric field in the two-stream in¬ 
stability in a short system compared to an infinite plasma 
can be very large. It is necessary to mention, however, 
that these simulations are carried out with certain simpli¬ 
fications similar to those made in the fluid model. In par¬ 
ticular, the sheath is not resolved, the ion background is 
immobile, and collisions with neutrals are omitted. The 
realistic sheath will allow some energetic plasma electrons 
to escape and may affect the structure of the wave inter¬ 
acting with the beam. If the ion dynamics is accounted 
for, the strong plasma oscillations may result in the mod¬ 
ulation instability which will create density cavities and 
affect the wave. Finally, the two stream instability can 
be suppressed by electron-neutral collisions if the colli¬ 
sion frequency is more than two times the collisionless 
growth rate. In very short systems, the temporal growth 
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FIG. 9. (a) Amplitude of the first maximum versus the rela¬ 
tive beam density in PIC simulations with L = 8.3A(, (green) 
and L = l.SAb (red). In (a), the solid blue curve is the theo¬ 
retical prediction (I22II for an infinite plasma, the dashed hori¬ 
zontal lines mark theoretical predictions given by Eq. (1^ for 
finite-length systems with L = S.SAb (green) and L = 1.5A(, 
(red), (b) Saturation amplitude versus the system length. In 
(b), the black curve is obtained with Eq. (1251) . the red vertical 
crosses and the green diagonal crosses mark first amplitude 
maxima in the PIC simulations with L = 8.3A(, (green) and 
L = l.SAb (red); the blue horizontal lines mark theoretical 
values (HU for the infinite plasma with the beam of relative 
density a = 0.00015 (dash), 0.0003 (short dash), and 0.0006 
(dash-dot). The PIC simulation values for L = 8.3A6 and 
L = l.hAb are obtained at points marked by vertical arrow 
B in Fig. [Ud) and by the vertical arrow in Fig. [Uh), respec¬ 
tively. 


rate is very small, which means that the neutrals present 
in a real beam-plasma system may simply prevent the 
instability from developing. 


VI. SUMMARY 

In summary, we have studied the development of the 
two-stream instability in a finite size plasma bounded by 
electrodes both analytically and making use of fluid and 


particle-in-cell simulations. We show that the instability 
reaches the asymptotic state when the wave structure has 
the same spatial profile and grows in time with a constant 
growth rate. The spatial structure of the wave is close to 
a standing wave but has a spatial growth along the beam 
propagation. We derived analytic expressions (I171I20I) for 
the frequency, wave number and the spatial and temporal 
growth rates. Obtained analytic solution agrees well with 
the values given by fluid and particle-in-cell simulations. 

The saturation of the instability occurs due to the over¬ 
taking of beam particles. Formulas for the estimate of the 
saturation amplitude (I241I25I) are derived and are in good 
agreement with the simulation results. The amplitude 
of saturation does not depend on the beam current but 
grows significantly for shorter systems. Compared to the 
value predicted for an infinite plasma, the saturation am¬ 
plitude for low-current plasma beam systems of length of 
a few resonance wavelengths may be higher by more than 
an order of magnitude. 
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